Virome analysis of Amblyomma americanum, Dermacentor variabilis, and Ixodes

24 A wide range of bacterial pathogens have been described in ticks, yet the 25 diversity of viruses in ticks is largely unexplored. In the United States, 26 Amblyomma americanum, Dermacentor variabilis, and Ixodes scapularis are 27 among the principal tick species associated with pathogen transmission. We 28 used high-throughput sequencing to characterize the viromes of these tick 29 species and identified the presence of Powassan virus and eight novel viruses. 30 These included the most divergent nairovirus described to date, two new clades 31 of tick-borne phleboviruses, a mononegavirus, and viruses with similarity to plant 32 and insect viruses. Our analysis revealed that ticks are reservoirs to a wide range 33 of viruses and suggests that discovery and characterization of tick-borne viruses 34 will have implications on viral taxonomy and may provide insight into tick35 transmitted diseases. 36 37 Importance 38 Ticks are implicated as vectors of a wide array of human and animal pathogens. 39 To better understand the extent of tick-borne diseases, it is crucial to uncover the 40 full range of microbial agents associated with ticks. Our current knowledge of the 41 diversity of tick-associated viruses is limited, in part due to the lack of 42 investigation of tick viromes. In this study we examined the virome of three tick 43 species from the United States. We found that ticks are hosts to highly divergent 44 viruses across several taxa, including ones previously associated with human 45 disease. Our data underscore the diversity of tick-associated viruses and 46 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom provides the foundation for further studies into viral etiology of tick-borne 47 diseases. 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom Introduction 70 Ticks (class Arachnida, subclass Acari) have been implicated as vectors in 71 a wide range of human and animal diseases worldwide (1-10). Approximately 72 900 species of ticks have been described and taxonomically classified into three 73 families: Argasidae (argasid or soft ticks) Ixodidae (ixodid, or hard ticks) and 74 Nuttalliellidae (11). Their propensity for feeding on on a wide array of hosts, 75 expansive range and long life cycle underscore the importance of tick 76 surveillance for the presence of potential pathogens. Argasid and ixodid ticks 77 combined transmit a greater diversity of viral, bacterial and protozoan pathogens 78 than any other arthropod vector (12). The worldwide incidence of tick-borne 79 disease is increasing, partly due to increased frequency of endemic tick-borne 80 diseases, as well as the discoveries of new tick-associated agents (13). 81 In the United States (US), bacterial agents are implicated in the majority of 82 tick-borne disorders. Lyme disease, caused by Borrelia burgdorferi represents 83 the most frequently reported tick-borne illness (14). Other bacterial agents, such 84 as Anaplasma, Ehrlichia, Rickettsia, other Borrelia species, as well as the 85 protozoan Babesia contribute to the overall spectrum of tick-borne disease (1, 6, 86 7, 14-16). Conversely, viral causes are diagnosed in only a fraction of tick-borne 87 disease cases (14). Despite considerable insights into the diversity of tick 88 bacteriomes, our understanding of tick-associated viruses is still limited. 89 Traditional viral isolation and identification methods using tissue culture have 90 isolated several tick-associated viruses but few have been characterized thus far. 91 In comparison to bacterial and protozoan agents, the literature associated with 92 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom tick-borne viruses in the Americas is also limited. Most of the literature to date on 93 suspected viral tick-borne pathogens focuses on those found in Europe, Asia and 94 Africa (3, 17). Powassan and Colorado tick-fever viruses have been historically 95 recognized as the only human tick-borne viral pathogens in the US (3, 18). 96 Recently, the Heartland virus, a new pathogenic tick-borne virus was 97 isolated from patients in Missouri and characterized by high-throughput 98 sequencing (HTS) (19). Heartland virus was shown to be phylogenetically similar 99 to severe fever with thrombocytopenia syndrome virus (SFTSV), a tick-borne 100 virus isolated from ticks and humans in China in 2009 (10). The emergence of 101 novel pathogenic tick-borne viruses as well as the dearth of data on tick viromes 102 suggests a need for viral surveillance and discovery in ticks. Although it seems 103 plausible that there are tick-borne viruses that have not been amenable to 104 isolation via tissue culture, to our knowledge, no extensive culture-independent 105 studies have been attempted to examine tick viromes. These studies may not 106 only identify viruses associated with acute disease but could also provide insights 107 into the pathogenesis of more controversial chronic illnesses associated with tick 108 bites (20, 21). Thus, to survey viral diversity, we examined viromes of three 109 human-biting ticks in the US (Amblyomma americanum, Dermacentor variabilis, 110 and Ixodes scapularis) by HTS. Our analysis and characterization identified eight 111 new viruses and indicates that ticks carry a wide array of previously 112 uncharacterized viral agents. 113 114 115 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom Materials and Methods 116 Adult I. scapularis, D. variabilis, and A. americanum were collected within a 2 117 square mile area of Heckscher State Park (Suffolk County, NY) in April 2013 118 (Figure 1). Ticks were pooled prior to nucleic acid extraction; two pools of I. 119 scapularis (N=30/pool), two pools of D. variablis (N=30/pool) and a single pool of 120 A. americanum (pool N=25) were prepared. A two-step purification and nuclease 121 treatment protocol was followed prior to extraction to enrich for viral sequences 122 and minimize host, bacterial and fungal template that can compete with virus 123 template in HTS (22). Each tick pool was homogenized in 500 μl of PBS, 124 followed by purification through .45 μM filter. The filtrate (237 μl) was treated with 125 1.5 μl RNase A for 15 minutes, followed by Turbo DNase (7.5 μl), Benzonase 126 (1.8 μl) and 2.7 μl of 1M MgCl2 for 45 minutes. All nuclease treatment was 127 performed at room temperature. 250 μl of the nuclease treated filtrate was added 128 to 750 μl of Nuclisens buffer and total nucleic acid (TNA) was extracted using the 129 EasyMag extraction platform (Biomerieux). TNA from each pool was eluted in a 130 35 μl volume. 131 132 Unbiased high-throughput sequencing 133 TNA (11 μl) from each tick pool was subjected to first and second strand cDNA 134 synthesis with Super Script III reverse transcriptase (Invitrogen) and Klenow 135 Fragment (New England Biolabs), respectively. Random primers (30ng/μl) 136 (Invitrogen) were used in both assays. Ion Shear Plus Reagents Kit (Life 137 Technologies) was used for double stranded cDNA fragmentation at 37°C for 25 138 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom min. Agencourt® AMPure® XP (Beckman Coulter) reagent (1.8 × sample 139 volume) was used for DNA purification. Ion Xpress Adapters and unique Ion 140 Xpress Barcodes (Life Technologies) were ligated to fragmented material by 141 using the Ion Plus Fragment Library kit (Life Technologies). Ligation was 142 performed at 25°C for 15 min and 72°C for 5 min. Ligated and nick repaired 143 products were purified with Agencourt® AMPure® XP (0.85 × sample volume) 144 and amplified according to manufacturers’ instructions with Platinum®PCR Super 145 Mix High Fidelity from the Ion Plus Fragment Library kit. Amplified products were 146 purified as previously described for the ligation reaction. Agilent® High Sensitivity 147 DNA Kit was used for library quantitation on the BioanalyzerTM 2100 instrument. 148 The concentration of each barcoded library was approximately 50nM. Libraries 149 were diluted to approximately 45nM and a pool of libraries in equimolar 150 concentrations was prepared. Ion OneTouch 200 Template Kit v2 (Life 151 Technologies) was used to bind barcoded libraries to Ion Sphere particles 152 (ISPS). Emulsion PCR of DNA linked ISPS was performed on the Ion 153 OneTouchTM 2 instrument (Life Technologies). Ion OneTouch ES instrument 154 was used to isolate template-positive ISPS. Ion PGMTM Sequencing 200 Kit v2 155 (Life Technologies) was used for sequencing of templated ISPS which were 156 loaded on the Ion 316TM Chip for further processing on the Ion Personal Genome 157 Machine® (PGMTM) System (Life Technologies). Approximately 600, 000 reads 158 were obtained for each library. 159 160 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom The de-multiplexed reads were preprocessed by trimming primers and adaptors, 161 length filtering, and masking of low complexity regions (WU-BLAST 2.0). The 162 remaining reads were subjected to homology search using BLASTn against a 163 database consisting of ribosomal and genomic metazoan sequences. Following 164 the processing, the remaining reads amounted to 116,946 (pool 1) and 131,676 165 reads (pool 2) for I. scapularis, 30,353 (pool 1) and 29,334 (pool 2) for D. 166 variabilis and 297,175 for the A. amblyomma pool. 167 168 The host-subtracted reads were assembled using the Newbler assembler (454, 169 v2.6). Contigs and singletons were then subjected to a homology search against 170 the entire GenBank database using BLASTn and the viral GenBank database 171 using BLASTx. Contigs and singletons with similarity to viral sequences from the 172 BLASTx analysis were again subjected to a homology search against entire 173 GenBank database to correct for biased e-values. For potential viral candidates, 174 close relatives were used to identify low homology regions in the genome from 175 BLASTx; gaps were filled in by PCR using primers specific to the assembled 176 sequence. Genome termini were obtained by 5’ and 3’ rapid amplification of 177 cDNA ends (RACE) (Clontech Laboratories). The final genome sequences were 178 verified by classical dideoxy sequencing using primers designed to generate 179 overlapping PCR products. 180 Tick screening 181 To determine the authenticity and prevalence of viral sequences in ticks, we used 182 cDNA generated from adult I. scapularis for a previous study as template for 183 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom PCR (23). The ticks were collected in 2008 from four separate locations in New 184 York State: two from Suffolk county (Heckscher State Park, Fire Island) and two 185 from Westchester County (Kitchawan Nature Preserve, Blue Mountain 186 Reservation) (Figure 1). cDNA from individual ticks was screened by PCR for 187 each novel virus identified by HTS. cDNA from individual virus-positive ticks were 188 used to generate viral genomic sequence equivalent to the sequences obtained 189 by HTS. The sequences were obtained by generating overlapping PCR products 190 using HTS sequences as a reference for primer design and sequence assembly. 191 All PCR products were verified by dideoxy sequencing. Sequences from 192 individual ticks were deposited in GenBank under accession numbers 193 KM048311-KM048322. All genome assemblies and alignments were performed 194 with Geneious v 6.1 and Mega 5.2 programs (24, 25). Phylogenetic trees were 195 constructed with Mega 5.2 using the maximum likelihood method with 1000 196 bootstrap replications. Amino acid trees were generated using the JTT matrix197 based model. Nucleotide trees were generated using the Jukes-Cantor model. 198 199 Results 200 Analysis of HTS data revealed the presence of eight previously uncharacterized 201 viruses (Table 1). The most diverse virome was observed in I. scapularis where 202 in addition to Powassan virus, we identified sequences representing six novel 203 viruses. Bunyavirus-like sequences represented more than 50% of the total 204 filtered reads in both I. scapularis pools, and were similar to viruses in the 205 Nairovirus and Phlebovirus genera. Two viruses were identified in each pool of 206 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom D. variablis; one virus was similar to phleboviruses and the other displayed 207 limited homology to insect viruses. Although the pool of A. americanum 208 generated the most reads of the five analyzed tick pools, a rhabdovirus was the 209 lone virus identified in this tick species. The complete genome characterization of 210 this virus, named Long Island tick rhabdovirus (LITRV) has been reported (26). 211 212 Arthropod genomes frequently contain integrated fragments of archaic RNA viral 213 genomes (27-29). The genome of I. scapularis contains numerous sequences of 214 viral origin, designated endogenous viral elements, or EVEs (27). To test for the 215 authenticity of the novel sequences obtained by HTS, we compared results of 216 RT-PCR vs DNA PCR as an additional indication of the template source i.e., 217 authentic viral nucleic acid vs. tick genomic DNA. Consistent with a model where 218 these novel virus sequences represent authentic RNA viruses and not EVEs, 219 products were obtained with RT-PCR but not with PCR. 220 221 Powassan virus 222 Two lineages of Powassan virus (POWV, family Flaviviridae, genus Flavivirus) 223 circulate in the US: lineage I, isolated mainly from I. cookei, and lineage II (also 224 known as deer tick virus), detected predominately in I. scapularis (30). Human 225 infection with viruses from either lineage has been linked with Powassan 226 encephalitis, a severe potentially life-threatening neurological illness (30, 31). 227 POWV lineage II has been detected in I. scapularis throughout the Northeast (23, 228 32, 33). In our HTS study, approximately 0.8% of all filtered reads from I. 229 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom scapularis pool 1 were from POWV lineage II. Assembly of all reads and contigs 230 from this pool resulted in assembly of 98% of the 10.8 Kb POWV genome. 231 Nucleotide (nt) comparison of the complete polyprotein coding sequence 232 (GenBank accession number KJ746872) indicated that the nt sequence of this 233 virus (designated strain LI-1) is 99.4% identical to isolate NSF001 (GenBank 234 accession number HM44059) obtained from Nantucket, MA in 1996 and 235 clustered with viruses in a subclade of POWV lineage II isolated from the 236 Northeast (Figure 2) (33). The nt identity to the subclade of lineage II made up of 237 isolates from the Midwestern United States was 93.4%. 238 239 Nairovirus 240 The genus Nairovirus (family Bunyaviridae) is comprised of 37 tick-borne viruses. 241 Its genome consists of three segments of negative sense single-stranded RNA, 242 designated small (S), medium (M), and large (L), that encode the nucleocapsid 243 protein (N), the envelope glycoproteins (Gn and Gc) and an RNA-dependent 244 RNA polymerase (L), respectively (34). Both I. scapularis pools contained 245 multiple contigs with sequence similar to viruses in the Nairovirus genus by 246 BLASTx. Assembly of these Nairovirus-like contigs, using CCHFV as a reference 247 genome, provided >90% coverage of the S and L segments. The complete 248 segments were obtained by overlapping PCR and 3’ and 5’ RACE. The 249 assembled sequences showed low similarity to other members of the genus, 250 suggesting that this virus, provisionally named South Bay virus (SBV) after its 251 geographic location, represents a novel nairovirus species. We also obtained the 252 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom complete S and L sequences from an individual SBV-positive tick (designated 253 SBV H38). Comparison of the genomic nt sequences obtained from HTS and 254 SBV H38 indicated they were 97% and 99% identical in the S and L segments, 255 respectively. SBV reads were the predominant viral reads obtained from both I 256 scapularis pools and accounted for 25% of filtered reads from pool 1 and 41% of 257 reads from pool 2. However, despite an exhaustive bioinformatics analysis, we 258 were unable to identify any contigs or reads with any similarity to Nairovirus M 259 segments. Attempts at virus isolation by inoculation of I. scapularis SBV PCR260 positive pools in Vero, Cos7, C6/36 and 297 cell lines were unsuccessful. 261 262 To ascertain whether the SBV sequences represented authentic viral reads, we 263 screened I. scapularis DNA and corresponding cDNA for the presence of SBV. 264 Partial S and L sequences were amplified in SBV-positive tick cDNA, but were 265 absent in the corresponding genomic DNA. Additionally, we found that I. 266 scapularis genome contains an integrated Nairovirus N gene open reading frame 267 (ORF) (accession number XM002414099), although it was not transcribed in any 268 tick sample. This sequence was more similar to SBV than other nairovirus 269 sequences, likely representing an ancestral integration of an SBV-like nairovirus 270 into the I. scapularis genome. 271 272 Conserved genus-specific sequences at the termini of each segment are a 273 feature of all bunyaviruses. In nairoviruses the termini sequence typically 274 consists of UCUCAAAAGA at the 5’ end. We found that the sequences of SBV L 275 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom and S termini were consistent with the Nairovirus sequence with the exception of 276 a single nt change at the fifth position in this sequence (UCUCUAAAGA). The 277 complete sequences of the L and S segments were deposited in GenBank under 278 accession numbers KJ746877-KJ746878 (originating from tick pools) and 279 KM048320-KM048321 (tick H38). 280 281 Phylogeny and L segment analysis 282 Of the 37 described viruses assigned to the Nairovirus genus, only seven are 283 completely sequenced and there is no published sequence data for 21 putative 284 members of the genus. For the remaining nine viruses, the only available 285 sequence data consists of a short (<450 nt) fragment within the L segment (35). 286 This region is part of the polymerase catalytic domain, and is highly conserved 287 with no nt insertions or deletions in any of the 16 nairoviruses analyzed to date. 288 Thus, we used the amino acid (aa) sequence from this region to determine the 289 phylogenetic relationships of SBV to the rest of the Nairovirus genus. SBV 290 contained a seven aa insertion in this region and had <45% identity to the next 291 closest virus, distinguishing SBV from other nairoviruses (Table 2). SBV did not 292 cluster with any of the described nairovirus serogroups, but fell outside of all 293 described viruses in this genus (Figure 3). 294 295 The L segment of SBV is 13892 nt long and contains a 13611 nt ORF that 296 encodes a 4536 aa protein. This represents the longest known nairovirus open 297 reading frame, over 600 amino acids longer than the L ORF of Dugbe virus. 298 on July 3, 2017 by gest http/jvi.asm .rg/ D ow nladed fom Comparison of SBV to L proteins from other nairoviruses revealed that SBV does 299 not cluster with any of the previously sequenced viruses, and contains only 28% 300 aa identity to the next closest virus (Erve virus) (Figure 4A, Table 3). In SBV, the 301 catalytic polymerase domain is located within the region between amino acid 302 (aa) residues 2600 and 3200 and contains all known viral RNA polymerase 303 motifs (pre-motif A, motifs A-E) (36). 304 In addition to the polymerase domain, the nairovirus L proteins may contain 305 several other protein motifs, such as ovarian tumor domain (OTU), 306 topoisomerase domain, zinc finger motif and a leucine zipper motif (36, 37). The 307 OTU domain represents a family of cysteine proteases and it is present in the N 308 terminus of all seven sequenced nairovirus L proteins. In CCHFV this domain 309 may function in modulation of interferon response in host cells (38). We found no 310 evidence of the presence of a functional OTU domain in SBV; the conserved 311 catalytic residues (D, C, and H in CCHFV) were absent and overall, the 312 SBV N-terminal region of the L showed little similarity to other nairoviruses. The 313 N terminal topoisomerase motif, zinc finger motif and leucine zipper motifs were 314